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Abstract 

Several groups have recently modeled evolutionary transitions from an ancestral allele to 
a beneficial allele separated by one or more intervening mutants. The beneficial allele can 
become fixed if a succession of intermediate mutants are fixed or alternatively if successive 
mutants arise while the previous intermediate mutant is still segregating. This latter 
process has been termed stochastic tunneling. Previous work has focused on the Moran 
model of population genetics. I use elementary methods of analyzing stochastic processes 
to derive the probability of tunneling in the limit of large population size for both Moran 
and Wright-Fisher populations. I also show how to efficiently obtain numerical results 
for finite populations. These results show that the probability of stochastic tunneling is 
twice as large under the Wright-Fisher model as it is under the Moran model. 

Keywords: Population Genetics, Stochastic Process, Stochastic Tunneling, Fixation 
Probability 



1. Introduction 



Evolutionary biologists have long understood that transitions between adaptive sets 
of traits may involve multiple substitutions separated by neutral or maladaptive inter- 
mediate states (Wright 1932 1. There has been a resurgence of interest in these ideas, in 



part because of advances in methods to measure epistatic interactions (e.g. Tong et. al 



2001 20041 and ability to observe evolutionary trajectories (Weinreich and Chao 20051. 
Several researchers have modeled evolutionary processes when epistatic interactions allow 
for multiple genotypes to have the same direct effect on fitness but experience different 



evolutionary dynamics because of differences in their genetic robustness) 


van Nimwegen 


et al. 1999 


de Visser et al. 


2003 


Proulx and Phillips 


2005 


Draghi et al. 


2010 


) or the 


local mutational landscape 


Wilke et al. 


2001 


O 'Fallon et al. 2007 


|. These scenarios 



can be called circum-neutral because alternative genotypes differ in their long-term evo- 
lutionary dynamics only because of the genomic circumstances in which they are found 



(Proulx and Adler 20101. 



Several groups have extended the theory to describe the rates and probability of 



transition along a multi-step evolutionary pathway. Weinreich and Chao ( 2005 ) took the 
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approach of calculating the total waiting time along various pathways and comparing the 



relative waiting times to reach a final state. Hermisson and Pennings (2005) considered a 



scenario where previously accumulated genetic variation may become adaptive following 
an environmental shift. In this scenario the population genetic dynamics of standing 
variation plays an important role in determining how evolution proceeds at the next step 
in the process (see also Kopp and Hermisson| (|2009)). Iwasa et al. (2004) derived approx- 
imate results on the waiting time and probability of a two-step sequence of mutational 



transitions using the Moran model, while Iwasa et al. ( 2003 ) derived results in a Wright 



Fisher model for a scenario where multiple mutations are required to escape the immune 
response. These results have been utilized by several other groups to study the rate of 



multi-step evolutionary processes (Durrett and Schmidt 



2008 


Lynch 


2010 


Lynch and 



Abegg 2010). Several other works have explored the probability and timing of multi- 



Weissman et al. 2009 Durrett et al 



step processes, as well as exploring the validity of approximations ( Schweinsberg 

' 2009] ) 



2008 



( 


2008 


) and 


Weissman 



et al. (2009) have presented branching process approximations for large populations that 
are equivalent to the large population size limit results for the Moran model presented 
here. 

The goal of this paper is to show how the finite population processes for both the 
Moran model and the Wright-Fisher model can be written and solved using the method 
of first step analysis. This helps to clarify some of the terms described by Iwasa et 
al. (2004) and gives an algorithm for efficiently solving the finite population Moran model. 
The Moran tunneling probabilities have previously been applied to Wright-Fisher pop- 
ulations without verifying that these results still hold. I show that the Wright-Fisher 
tunneling probabilities differ from the Moran probabilities by a factor of 2. This correc- 
tion will allow stochastic tunneling results to be applied to a wider range of scenarios. 
I also compare the large population size approximations for the rate of tunneling with 
simulations and exact calculations for small population size. 

1.1. Preliminary definitions and results 

By considering the population level evolutionary process as a series of transitions 
between populations fixed for a single genotype we can calculate the waiting time for 
the population to become fixed for secondary mutations. So long as Nfj, << 1 we 
will seldom have multiple mutants arising in the same generation. This approach also 
assumes that each attempt at tunneling, if unsuccessful, is over before another primary 
mutant arises. Determining when this conditions actually holds is more difficult because 
the sojourn time of the primary mutant goes up as its selective disadvantage decreases. 
In the case of circum-neutral primary mutants, the sojourn times are characterized by 
large variances that become undefined as population size approaches infinity. A rigorous 
analysis of the parameter combinations that allow this approximation to be applied is 
provided in |Schweinsberg (2008) and Weissman et al. ( |2009 ). 

The first mutational step (the primary mutant) is assumed to have relative fitness r < 
1, while the second mutational step (secondary mutant) is assumed to have fitness a > 1 
relative to the ancestral allele. In the case where r is exactly one, the first mutational 
step has no direct effect on fitness and the primary mutants can be considered circum- 



neutral (Proulx and Adler 2010 1. Such circum-neutral substitutions do not directly affect 
reproductive fitness but do alter the long-term evolutionary trajectory of the population. 
The ancestral population can evolve to be fixed for the secondary mutant either through 

2 



a sequential mutational pathway or because a lineage of primary mutants destined for 
extinction produces a secondary mutant which is destined for fixation, a process termed 



stochastic tunneling by ( Komarova et al. 2003 1 . 



The waiting time until a secondary mutation becomes fixed can be expressed in terms 
of the waiting times for the sequential and tunneling paths. I define the per generation 
probability of successful sequential substitutions Si and 5*2 and the per generation prob- 
ability of the opening of a successful tunnel as T . The waiting time for the transition 
between population states is well described by an exponential waiting time so long as 



population size is not too small (Iwasa et al. 2005). This means that the process is 



characterized by a race between waiting for a primary mutation to arise and fix and the 
start of a tunneling pathway. The expectation of the total waiting time until a secondary 
mutation is given by 

FM _ T Si(Si+S 2 + T) 

lJ (T + Si) 2 S 2 {T + Si) 2 ' K) 
where the first term represents the contribution to the expected waiting time from tunnel- 
ing pathways and the second term represents the contribution from sequential pathways. 
If T = this is simply the sum of the waiting times for primary and secondary mutations 
to sequentially fix. This approximation ignores the time that it takes for beneficial mu- 
tations to spread through the population and the amount of time that primary mutants 
are segregating before a secondary mutation arises. The time required for alleles destined 
to fix to spread to fixation is typically much smaller than the waiting times for them to 



arise, and in any case it can be simply added to the total waiting time (see Lynch and 



Abegg 2010). 



The per generation probabilities of sequential fixation are 

51 = NtnU(r), (2) 

5 2 = Nfi 2 U(a/r), (3) 

where N is the haploid population size(for simplicity I assume this is approximately the 
effective population size as well), fii is the probability that an ancestral allele will mutate 
into a primary mutant, is the probability that the primary mutant will mutate into 
a secondary mutant, and U(x) is the fixation probability of a mutation with relative 
fitness x when initially present as a single copy. Because this follows sequential fixation 
of mutants, the secondary mutant is invading into a population fixed for the primary 
mutant, giving it a relative fitness of a/r. 

Following Iwasa et al. (2004) , the probability of tunneling can be written as 

T = (1 — U(r))(l — E[no secondary substitution|extinction]), (4) 

where E[no secondary substitution|extinction] represents the probability that no success- 
ful secondary mutations arise while the primary mutant is segregating conditioned on 
the eventual extinction of the lineage of primary mutants. This can be related to the 
unconditional expectation by 

£7 [no secondary substitution! extinction] = E[no secondary substitution] (1 — U(r)) (5) 



(Iwasa et al. 2004 ). This provides a simple relationship between calculations made using 
the conditioned trajectory of primary mutations and the unconditioned trajectory of 
primary mutations. 
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2. Moran Model 



The Moran model (Moran 1962) follows a population of size N described by a vector 
x where Xi indicates the number of individuals of genotype i. In the Moran process, each 
unit of time either 2 elements of x change by one unit each in opposite directions or x 
remains constant. This model is often conceptually described as one where an individual 
is chosen to reproduce at random weighted by their reproductive output. Population size 
is kept constant by choosing one of the original population members to die (it may be the 
one that reproduced). Mutation causes offspring to differ from their parent's genotype 
with a probability defined by the mutation rate. The scale of time in this model is 
population size dependent; a generation is measured in terms of N time steps. 

Because the number of individuals of each genotype can change by at most one unit 
each time step, the Moran model can be expressed as a Markov chain whose matrix 
definition is tridiagonal. This is in sharp contrast to the Wright-Fisher model where 
any population state can move to any other state in one generation (albeit with low 
probabilities) . This simple matrix structure allows many features of the stochastic process 
to be expressed algebraically. 

The evolution of the population can be described as a series of transitions between 
states described by the complement of segregating alleles. Assuming that primary mu- 
tations occur rarely enough so that multiple primary mutants do not typically arise 
together, the transition probabilities between states can be based on the introduction of 
a single individual mutant. As this approximation breaks down more error will be in- 
troduced into the transition probabilities. The ancestral population is monomorphic for 
the ancestral genotype. Following the introduction of a primary mutant, the population 
will evolve with two genotypes for some time until either the lineage of primary mutants 
goes extinct or a secondary mutant lineage arises and does not go extinct. 

Following the introduction of a primary mutant, the population is composed of i 
primary mutants and N — i ancestral alleles. In the absence of mutation, the Markov 
transition probabilities are 

Pr(i — 

Pr(i -> i) 

Pr(i + 



i(N-i) 


N(ir + (N - 


-i)) 


i 2 r + (N - 


^f 


N(ir + (N - 


-i)) 


ir(N — i 


) 



N(ir + (N-i))' 

where < i < N and r is the relative fitness of the primary mutant. Note that this 
model ignores the change in the number of primary mutants due to their mutation into 
secondary mutants. Following the introduction of the primary mutant, it may produce a 
lineage that eventually goes extinct or eventually becomes fixed in the population. It is 
useful to describe the population process conditioned on the eventual extinction of the 
primary mutation in order to consider these two scenarios separately. The conditional 
process can simply be described by 

Pr(i — > j|extinction) = Pr(i — > j) — 
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where tt* represents the the probability that the lineage goes extinct given that there are 
currently i mutants in the population ( Ewens 



I use first step analysis (Taylor and Karlin 



1973) 



1984) to in order to find the total proba- 



bility that no successful secondary mutant is spawned from a lineage of primary mutants 
destined to eventual extinction. Let Vi be the probability that no successful secondary 
mutants are spawned from a lineage beginning with i mutants. Vi can be implicitly 
defined as the probability that no successful mutants arise in the current time step mul- 
tiplied by the probability that no successful mutants are produced in the future. For the 
process conditioned on eventual extinction of the primary mutant lineage we have 

Vi = (1 — cor — ) ( Pr(i -> i — 1) % ~ l Vi-i + Pr(i — ¥ i)—Vi + Pr(i -> i + 1) %+1 Vi + \ j , (6) 

N \ TTi Hi TTi J 

where the composite parameter to = fi^Uia). Note that vq = 1 and ttn = 0. This then is 
a system of N — 1 linear equations with N — 1 unknowns. The probability of tunneling 
is simply T = Nfj,i(l — i>i)7Ti. 

2.1. Algorithm for solving the finite population size model 

For finite populations the system of equations can be represented as a tridiagonal 
matrix. This system can be solved numerically using a mathematical computing pack- 
age. Many results are known for the matrix inverse and eigenvectors of tridiagonal 
matrices (Usmani 1994 da Fonseca 2007). These results can be used to numerically 



calculate the eigenvectors for even large population size because they involve less than 
AN multiplications and additions. Because they use recursive calculations there is no 
constraint imposed by memory levels and computation time basically increases linearly. 
On a MacPro with a 3.2 GHz Xeon processor running Mathematica 7 the calculation for 
population size of 10 6 takes about 50 seconds. 

Given our system of equations ([6| we can write a matrix equation of the form 

Av = x, 



where v is the vector of Wj and 

/ax &i 



ci a 2 b 2 
C2 a 3 

C3 



V 





on -2 

CjV-2 CtJV-l/ 



The elements of the matrix can be found from equation (|6| and I provide formulas for a^, 
bi, Ci and x\ below. Note that this system has N — 1 rows because it is conditioned on 
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the eventual extinction of the primary mutant, (da Fonseca 2007) defines the recursive 
equations 













= 1, 


01 


= »i 








= 1, 


N-l 


= G-./V— 1 ■ 



he, 



Each element of A -1 can be expressed as an algebraic expression of 8 and </>. Because V\ 
is the first element of A _1 x and because only the first element of x is non-zero we need 
only calculate Aj"J. 



A^ = |^. (7) 



Given the system of equations ^ we have 



h = 



l 



{N - 1)(N -ujt) 
N 2 (N - l + r)(l-vri) 
(N - if + i 2 r N — iru 
N(N-i + ir) N 

(N — i)ir(l — Tti+i) N — iruj 
N(N -i + ir)(l-ni) N 
_ (N — i — l)(i + 1)(1 — 7Ti) iV - (i + l)rw 
Q ~~ N(N -i + ir)(l-n i+1 ) N 

Finally, the total probability that no successful secondary mutations are produced while 
a lineage descending from a single primary mutant is extant is 

= _A_ (_ f-i _ ™\ N - 1 \ ( H ) 

Vl 6 N -i\ \ N I N(N — 1 + r)(l — 7Ti) / lj 

This method can be used to numerically solve for v\ and is reasonably quick even in large 
populations. 

3. Wright-Fisher Populations 

3.1. Finite Population Size 

In the Wright-Fisher formulation the population at generation t + 1 is found by 
sampling gametes produced in generation t. So long as the number of gametes produced 
is reasonably large, the probability distribution for adults in generation t + 1 is binomial 



G 



such that 



Pr(0 -> 0) = 1, 

Pr(i -> j) = 
Pr(A ->• JV) = 1, 



N 



j ) \%T + (N - i) 



1 - 



ir + (N - i) 



N-j 



represents the probability that the population goes from i mutants to j mutants in 
one generation. Again I use a first step analysis to calculate the probability that no 
secondary mutations arise beginning from a single mutant. Exact calculations of the 
fixation probabilities for the finite Wright-Fisher model are not available, so this system 
cannot be converted to the process conditioned on eventual extinction of the primary 
mutant. This means that the probability of tunneling will have to be back-calculated 
from equation (JSJ. If the primary mutant becomes fixed then the probability that a 
successful secondary mutation is spawned is 1. 

I define v, as the probability that no successful secondary mutations are spawned 
starting from the state where i primary mutants are present. For each possible state 
i, the probability that no successful secondary mutants are produced is the probability 
that none of the i primary mutants immediately produce a successful secondary mutant 
((1 — uj) 1 ) multiplied by the sum of the probabilities that the next generation contains 
j primary mutants multiplied by the probability that a lineage starting with j mutants 
never produces a successful secondary mutant. This is slightly different from the Moran 
model where only one secondary mutant can possibly arise at each time point. This gives 



vo = 1, 



N 



Vi = (1 -u) l ^Pr{i -> j)vj, 

3=0 

v N = 0. 

The equations for 1 < i < N — 1 can be rewritten as a sum of Vj terms as follows 

N 

= (l-w) i ^Pr(i^j)u j -^ 



N 



= (i - u>y ^2 Pr ^ ~^ + Pr (* _ l ) + ^2 Pr ^ ~^ fi v 3 

\3=0 j=t+l 

This is a system of linear equations in Vj and can be written in matrix form as Av = x 
where 

/ 1 ... ... 0\ 

(1 - w) 1 Pr(l 0) (1 -u^PKl 1) - 1 (1 -w) 1 Pr(l 2) 

(1 - u) 2 Pr(2 -> 0) (1 - u) 2 Pr{2 -> 1) (1 - uj) 2 Pr{2 -> 2) - 1 



V 



1/ 



V2 



\v~nJ 



W 



The solution to v = A _1 x can be found numerically using standard mathematical pack- 
ages. Compared to the solutions for the Moran model, the Wright-Fisher model requires 
many more computational operations for the same population size. 

Recall that v\ is the unconditioned probability that no successful secondary mutations 
are produced from a lineage of primary mutants descending from a single primary mutant. 
We can calculate the probability that no successful secondary mutations are produced 
conditioned on the eventual extinction of a lineage of primary mutants descending from a 
single initial mutant using the approximate fixation probability for a single copy mutant 
(Crow and Kirnura 19701 and equation |5]) as 



vi = vi 1 - 



1 - er 2{r -^ 

1 - e -2W(r-i) 



4. Large Population Size Approximations 

In very large populations, the dynamics of newly introduced mutants can be modeled 
as a branching process. In this limit, segregating mutants do not interact and mean fitness 
is not altered by their spread in the population. The probability that a newly introduced 
primary mutant leads to the eventual maintenance of a secondary mutant can then be 
calculated based on this branching process, while fixation probabilities of beneficial alleles 
may still be modeled using finite population size results. In other words, population size 
enters the calculations in two ways; the finite population size causes frequency dependent 
interactions between segregating mutant alleles and causes the fixation probability of 
beneficial mutants to increase. The large population size approximation takes this first 
population size to be large while leaving this second population size finite. 

The main feature of the branching process approximation that allows this to be 
calculated is that the probability that no secondary mutations persist as a function of 
the number of initial mutants scales as Vi = a 1 , where a represents the probability that 
no secondary mutations persist when a single primary mutant is introduced. Once a has 
been solved for the per generation probability of tunneling is described by 

T = N f x 1 (l-a)(l-U{r)). (9) 



4.1. Moran Model 

The calculation of the expected probability of a successful secondary mutation in the 
large population limit begins with the same first step analysis equations as in a finite 
population. Starting with equation ^ setting i = 1 and taking the limit as N — > 00 we 
have 

V2 = + (1Q) 
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Inserting = a z M into equation ( 10 ) gives 



2r 



(11) 



Note that this is slightly different from the formula Iwasa et al. (2004) give be- 
cause theirs is an expression in terms of the unconditional process. Recalling that 
E[P | extinction] = B[P](1 — U(r)) their expression can be rewritten in terms of the 
conditional process as 



aiMN 



2- v /(l-r) 2 +2(l + r)7 
1 + r 



(l-f/(r)). 



(12) 



In general, these equations are extremely similar and will often be so close as to be 
numerically indistinguishable for any empirical applications. Notable exceptions are 
when the population is small enough that the Iwasa et al. (2004) breaks down and 
«imn > 1- 

In the case of a circum-neutral primary mutation (r = 1) the two expressions become 



a M 
«imn 



1+ 2 



- (u - V^(4 + ^)) 



(1 



N — 1 
N ' 



Recalling that ui = fj>2 ^_ ^j^ N and taking the limit as N — > oo we get 



lirn/v_ 
limjv^. 



i «M 



«imn 



(1 



a — 1 NO -, a — 1 

-M2) 2 - 1 



2a 



2a 



/'2 



a- 1 



-1*2 



(13) 
(14) 

(15) 
(16) 



(See ( Schweinsberg 20081 for an alternate derivation). These expressions are most dif- 
ferent when a is large and fi2 is large. However, even for an unrealistically large value 
fj,2 = 10 -2 the expressions are different by less than 1%. 



4.2. Wright-Fisher Model 

Using the branching process approximation we have Vj 
this gives 



a 



WF- 



For 1 < i < N - 1 



/v 



(i-«)'E 



N 



j J \ir + (N - i) 



1 



ir + (N - i) 



N-j 



"WF- 



(17) 



For any specific i, Vi can be approximated by noting that the binomial probabilities 



Ross 


1988 


Iwasa et al. 


2003) 



equation (17 1 in the limit of large N using the Poisson approximation gives 

= (l-w)* (e-Ki-owpjy 



v i = (l-u>) il £ 

j=0 



"(ir) j 



"WF 



(18) 



A similar result was obtained by Iwasa et al. (2003) for multi-step pathways involving 
multiple routes to a higher fitness mutant. For a primary mutant lineage starting with 
1 individual the probability of tunneling is the solution of 



a WF = e- r < 1 - awF '(l- W ). (19) 
This implicit definition of «wf can be easily solved numerically for specific values of r 



and w (figure B . 1 ) 



The implicit equation for the probability that no successful secondary mutants are 
produced has a convenient interpretation. In large Wright-Fisher populations the dis- 
tribution of offspring produced by a single parent is Poisson. A newly arising primary 
mutant has probability oj of immediately mutating and spawning a lineage of secondary 
mutants that is destined to fix. If it docs not (with probability 1 — oS) then it produces 
a Poisson distributed number of offspring mutants with mean r, each of which has an 
effectively independent probability 1 — %p of spawning a secondary mutant lineage itself 
destined to fix. The probability that none of these primary mutant offspring spawn a 
successful secondary mutant lineage is simply the term from the Poisson distribution 
with parameter r(l — owf)- Thus, awF is equal to the product of the probability that the 
lone primary mutant does not immediately produce a successful secondary mutant with 
the probability that none of the primary mutant offspring produce a successful secondary 
mutant lineage. 

4-3. Comparisons between the models 

Define Tm and Twf as the probability that a successful secondary mutation arises 
from a lineage of primary mutants founded by a single primary mutant in the Moran and 
Wright-Fisher models, respectively. 

For the Moran model where r = 1, Tm = \/w(l + w/4) — % (i.e. 1 — ctyi from 



equation (13)). For small to this can be approximated using a series expansion (see 
appendix Appendix A ) . Likewise the Wright-Fisher model can be approximated using 



equation (19) and again expanding around small uj (see appendix Appendix A) 



t m «^-| + om 3/2 , (20) 

T WF «y2^-^+OM 3 / 2 . (21) 

Noting that the fixation probability approaches U(a) = (a— 1) in large Moran populations 
and U(a) = 2(a — 1) in large Wright-Fisher populations and substituting u> = 2/X2?7(a) 
gives 

Tm^^^-I)- ^"^ , (22) 
4/i 2 (a- 1) 



r WF ^2 v /M2(a-l)- ^ 2l 3 ; - (23) 

These calculations show that for the same level of \Xi and a the Wright-Fisher tunneling 
probability is a factor of 2 larger than for the Moran model. 

In both models the probability of tunneling depends on the distribution of the num- 
ber of primary mutants spawned by a lineage founded by one primary mutant and the 
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curvature of the function relating the probability a successful secondary mutant arises to 
the total number of primary mutants in a given lineage. In the limit of large population 
size, where a branching process approximation is valid, this distribution can be exactly 
determined (Dwass 1969). Examining the relationship between the expectation of the 



probability no successful secondary mutant arises and the distribution of the primary 
mutant lineage size can help to understand this effect. Once population size gets to be 
reasonably large, the distribution of primary mutant lineage size remains almost constant 
except at the tail. In the Wright-Fisher case, when r — 1 we have a critical branching 
process. Conditioned on eventual extinction of the primary mutant lineage, we can calcu- 
late the expected number of primary mutants spawned in a lineage arising from a single 



primary mutant (see Appendix B ) . Alternatively, we could calculate the distribution of 
primary mutant lineage sizes and then calculate the total probability that no successful 
secondary mutations are spawned. Figure |B.3| shows how the curvature of the function 
for the probability of successful secondary mutations as a function of primary mutation 
lineage size would lead to a reduced total probability that a secondary mutation becomes 
fixed. Because of Jensen's inequality, the expectation of the function of the random vari- 
able is larger than the function of the expectation. Because the probability cannot be 
below 0, no matter how large the primary mutant lineage, the contribution of very large 
primary mutation lineages is negligible. Thus, T is considerably smaller than we would 
predict based only on the expected number of primary mutants spawned by a lineage 
destined to eventually become extinct. 

For cases where r < 1 approximating around small uj is more straightforward (see 
appendix Appendix A ) . Approximating and substituting in the values for w gives 



T M « /j,(a - 1) 
T WF « 2fi(a — 1 



r 



1 - r 



(24) 
(25) 



For r near one this approximation fails breaks down and either the exact equation or the 
neutral approximation should be used (figure B.l). Again, the probability of tunneling 



is larger for Wright-Fisher populations, but now they are different by a factor of 2r. 
The factor of 2 is again because of higher fixation probabilities in the Wright Fisher 
model while the factor of r is due to the Moran model assumption that mutations only 
occur during reproduction. If, alternatively, mutations were assumed to happen at a 
constant rate per generation (where generations encompass N elementary steps of the 
Moran process) , then the results would only differ because of their fixation probabilities 
(to the first order approximation). 

These results agree with the Moran population results of Iwasa et al. (2004). Iwasa 
et al. (2004) characterize the regime where r < 1 by arguing that it arises as the prod- 
uct of the equilibrium frequency of the deleterious primary mutation, the probability 
a secondary mutant arises, and the probability of successful fixation of the secondary 
mutation. The results presented here provide a different interpretation. The probabil- 
ity of tunneling as derived here is the average over independent trajectories following 
the introduction of a single mutant. Such trajectories are never at equilibrium, but do 
spawn a characteristic distribution of mutants before their lineage goes extinct. In both 
Moran and Wright-Fisher populations, deleterious mutants produce an average of 
descendants (see Appendix B ) . An important caveat is that equations ( 24 ) and ( 25 1 only 
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apply when r < 1 and w is small. In these cases equations ( 11 ) and ( 19 ) should be used 
instead. 



5. Multiple intermediate steps 



This approach readily extends to the case where multiple intermediate mutational 
steps stand between the ancestral state and any mutants that have improved fitness. 
Such multi-step probabilities have been calculated for the Moran model by Iwasa et 



(2010 


); 


Lynch 


(2010 


), and 


Weissman 



proximation is valid when /i is small relative to the inverse of the population size squared. 



Lynch (2010) considered both the multistep tunneling probability for both neutral and 



deleterious intermediates but did so by assuming that each intermediate deleterious mu- 



tation first rose to mutation selection balance. Additionally, Iwasa et al. (20041 consid- 
ered a immune response escape mutants in a multi-step Wright-Fisher model. Weissman 
et al. ( 2009 1 examined the probability of multi-step tunneling for an arbitrary number 



of intermediate mutations and rigorously defined the regions of parameter space where 
the stochastic tunneling approximation is valid. They note that the stochastic tunneling 
approximation, as used here, is only valid when population size is large. Specifically, the 
required population size goes up if there are many intermediate mutations and they are 



only weakly selected against ( Weissman et al. 2009 1 . 



The probability that a multistep tunnel is opened can be calculated recursively. Using 



equation (11 1, but writing u> generically as fiU, I define a function 

1 + rfaU + 1) - v/(l - r) 2 + (2 + HJJXJ + 2))r/iU 



g(p,U) = l- 



2r 



(26) 



For a 2-step process, the probability of a tunnel is simply g(fi2, U(a)). Further recursions 
give the total probability for longer tunnels, so that a three-step tunnel opens with 
probability g(p 2 , 9 (^3, U(a))). 

In the case where the sequence of mutations all have the same mutation rate we 
think of the probability a fc-step tunnel is opened by suppressing the variable /j, in g and 
recursively applying g to U (a) . The solution can be found graphically by cobwebbing 



the graph of g (figure B.2| (Adler 1998). Longer tunnels have lower probabilities, and 
the decrease in probability depends on the shape of g. If g'(0) < 1 then probability of 
more complex tunnels opening decreases towards (because g has no fixed point), but 
if g'(0) > 1 then the adding more intermediate mutations has a decreasing effect on the 
tunneling probability. Note that g'(0) < 1 if r < jj—. 

For Wright-Fisher populations the probability of tunneling was described implicitly 
as a solution to a transcendental equation. However, the same recursive approach can 
be taken to find successive tunneling probabilities. Substituting u> = fiU into equation 
( 19 ) and noting that at a fixed point of the recursion a = 1 — U gives 

1-U = e- r(c/) (l- fiU), (27) 

where solving for U gives the fixed point of the recursion. Decreasing r decreases the 
fixed point. Approximating around U = and solving gives 

2(/i + r - 1) 



r(2/x. 
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•r) 



(28) 



where Uoo represents the fixed point of the recursion. If r < 1 — fx then there is no positive 
fixed point. This fixed point represents the infinite recursion for the probabilities and 
relies on the assumption that the time-scales of fixation and mutation can be treated 
separately. However, Weissman et al. (2009) found that these conditions narrow as the 
the length of the pathway considered increases. As the length of the pathway increases, 
the probability of a series of sequential fixation events increases. However, when equation 
(28) has no fixed point we can still infer that tunneling across long valleys will not occur. 



6. Simulations of finite population 

The approach taken by Iwasa et al. (2004) involved first approximating the Moran 
process by a small time-step approximation and then using special functions and heuristic 
arguments to arrive at an approximation for the rate of tunneling. This method implic- 
itly assumes large population size and explicitly ignores some higher-order terms. My 
approximation explicitly assumes large population size to derive a result for the limit as 
population size goes to infinity. Numerical solutions can be used to evaluate the ability 
of these large-population size approximations to predict the rate of tunneling in finite 
populations. I simulated the discrete time Moran and Wright-Fisher models in order to 
assess the accuracy of the different approximations. My simulation draws waiting times 
for mutations and then tracks individuals in populations while multiple mutations are 
segregating. Once the secondary mutation has reached a significant size its fixation is 
virtually guaranteed, and the simulation is stopped. 

For the Moran model, the Iwasa et al. (2004) solution and my approximation are 



displaced from the numerical solution in opposite directions (figure B.4). The Iwasa et 
al. (2004) approximation overestimates the waiting time to a secondary mutation. Both 
approximations do extremely well once population size is larger than about 100 (for the 



parameter values in figure B.4). For smaller population sizes the variance in the wait- 
ing time is so large that, in practice, it would be hard to distinguish the alternative 
approximations. It is interesting to note that, in terms of displacement from the ac- 
tual waiting time, the Iwasa et al. (2004) approximation performs better, even though 
it intentionally ignores some terms. This is apparently because the large-population ap- 
proximation underestimates the waiting time (because mean fitness is not altered by the 
spread of mutants in the large-population limit), and by chance the terms that Iwasa et 
al. (2004) exclude happen to push the approximation further off. 

For the Wright-Fisher model, I compared my solution with the one presented by 



Lynch and Abegg [Lynch and Abegg ( 2010| ) (figure B.5). Lynch and Abegg applied the 



Iwasa et al. (2004) approximation to Wright-Fisher populations but modified the calcu- 
lations to adjust for the possibility of multiple mutational hits in a single generation. My 
approximation for the rate of tunneling in Wright-Fisher populations includes a factor 
of \/2 that is left out if one simply applies the Iwasa et al. (2004) Moran approximation 
to Wright-Fisher populations. At small population sizes, the sequential fixation pathway 
dominates and the two approaches yield similar predictions. At intermediate population 
size the predictions are most different, but still have qualitatively similar patterns. For 
the parameters used in figure |B.5[ the two predictions always differ by less than about 
0.3% and are within about 8 million generations of each other. For intermediate popula- 
tion sizes, the finite population matrix solution does capture the behavior of the system. 
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As population size grows larger the simulated results move towards my approximation, 
just as the matrix solution converges towards my approximation. 



7. Conclusions 

Multi-step evolution is becoming more widely recognized as an important component 



of the evolutionary process ( 


Weinreich and Chao| 


2005 


|Hermisson and Pennings 2005 


Kopp and Hermisson[ |2009| | 


Durrett and Schmidt) 


2008 


|Lynch[ 2010| Lynch and Abegg 



(Iwasa et al. 



berg 2008 Weissman et al. 2009 1 and for specific instances of the Wright-Fisher model 



2003 1 . The Moran model results have been used in models of Wright-Fisher 



populations (Lynch and Abcgg 2010) and the distinction between the two models has not 
received much attention. The results I present here are derived using elementary meth- 
ods from stochastic processes theory. I present exact expressions for the limiting case of 
large population size in both Moran and Wright-Fisher populations. Some methods for 
efficiently numerically evaluating the finite population size solutions are presented. 

The Wright-Fisher and Moran models differ in two different but related ways. First, 
the branching process calculation of the probability of tunneling depends on a sum over 



the number of mutant offspring produced by a single mutant (compare equations ( 11 ) and 



(17)). For the same mean selection against primary mutants, the variance in the number 
of offspring under the Wright-Fisher model is 1/2 as large as under the Moran model. 
The tunneling probability also depends on to, the product of the secondary mutation 
rate and the probability of fixation of the secondary mutation. Second, the fixation 
probability also depends on the distribution of the number of offspring and again differs 
by a factor of two between the models. For the same increase in fitness, the probability 
of fixation is twice as large under the Wright-Fisher model as under the Moran model. 
Because these are both introduced inside a square-root function, the total effect is that 
the probability of tunneling is twice as large under the Wright-Fisher model as compared 
to the Moran model. The general prediction is that when the offspring distribution has 
lower variance then the probability of tunneling will increase. 

For Wright-Fisher populations, the probability of tunneling is the solution to an 
exponential equation which has a simple intuitive explanation. The total probability 
that no successful secondary mutants are produced is the product of the probabilities 
that the initial primary mutant does not immediately produce a successful secondary 
mutant, 1 — u) and the probability that none of its primary mutant progeny produce a 
successful secondary mutant. 

My analysis shows that tunneling in the Wright-Fisher model is more likely than in 
the Moran model. Stochastic simulations show that the Wright-Fisher approximation 
does indeed capture the mean behavior of the evolutionary process once population 
size is relatively large. The improvement over applying the Moran approximation to 
the Wright-Fisher scenario is quite minor, but there is no added difficulty in using this 
correct approximation in the future. 
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Appendix A. Approximations around small oj 

The probability of tunneling in a large Moran population is 

T M = v/w(l+w/4) 



ui 
2 ' 



(A.l) 



The second term {—%) is already linear and does not need further approximation. We 
would like to approximate this function for small a; as a power series of terms aW 2 . Define 
/(£>) = (1 + Cj) and write the radical as y/u}f(u)), where Cj will later be set equal to to. 
We will construct a Taylor's series approximation around Cj = 0. Note that /(0) = 1 and 
/'(0) = 1/4, and all higher derivatives of / are 0. The series is as follows 



\Juf(Cj) sa ^/uj + Cj 



c/'(0)(w/(0))- 1/2 0? (wf (0)) 2 ( W /(0))- 3 / 2 



2! 



2 2 



(A.2) 



Taking just the zero order terms gives 

r M «Vw-|+0[o>]. (A.3) 
After setting cD = w, the general expression for the ith term in the Taylor's expansion 



is 



L 1 ! — I (W(i/2^0. 



(A.4) 



This series converges so that s * = V^\/^ + w/4. So a good approximation is 

T M «\Wl W 4 -f- ( A -5) 

For the Wright-Fisher model, when r = 1, a is represented implicitly by a = e 1_a (l — 
cj). Recall that Twf = 1 — a and express it in terms of Lambert's W (the solution of 
W{z)e w ^ = z) gives 

T WF = l + W{^-). (A.6) 



Many methods of approximating W are known, we use the power series expansion pre- 

(A.7) 



sented by Corless et al. (1996), 



W{z) « -l+p-l/V + 



where p = y / 2(ez + 1). Substituting this back into equation (A.6) gives 



Twf 



2uj- -uj + 0[uj} 1/2 . 
3 



(A. 



For situations where r < 1 a similar approach can be used. From equation (JllJ) we 
have 



T- 



M 



1 + r(w + 1) - v/(l-r) 2 + (2 + r(w + 2))7 
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2r 



(A.9) 



Away from r = 1 this can be approximated using a Taylor's series expansion to give 

Tm^-^+OM 3 - (A.10) 
1 — r (1 — r) 6 

For the Wright-Fisher model with r < 1 a regular perturbation approach can be used 
to the probability of a tunnel opening. First write a as a function of u giving 

a((j) = e- r(1 - Q(w)) (l-a;), (A.ll) 

and note that a(0) = 1. Differentiating the equation with respect to u> gives 

a'{w) = ra'(o;)e- r(1 ~ Q( " )) (l - w) - e - r (i-«H). (A.12) 

Setting u> = gives 

a'(0) = ra'fa;) - 1 a'(0) = — . (A.13) 

1 — r 

Thus, 

Twf = 1 — a:(w) ~ ^ — = 2u(a — 1) — - — . (A.14) 
1 — r ' 1 — r 

Appendix B. Mean number of primary mutants spawned by a lineage des- 
tined to become extinct 

Appendix B.l. Moran populations 

First step analysis can be used to calculate the mean number of mutants spawned 
by a lineage descending from a single initial mutant. Again I condition on the eventual 
extinction of the mutant linage. Note that in the Moran model, time is scaled by the 
population size and this must be kept in mind when using results based on the number 
of mutants produced to estimate the probability of a secondary mutation. The first step 
analysis yields the system of equations 

A = Pr(l -> 0) — (1 + D ) + Pr(l -> 1)(1 + D x ) + Pr(l -> 2) — (1 + D 2 ) (B.l) 

7Ti 7T1 

A = Pr(i -> i - l)^ 1 ^ + A-i) + Pr(i -> + A) + Pr(i -> i + l)^ 1 ^ + A+i 

(B.2) 

Av-i = Pr(iV — 1 — > A — 2)^^(A - 1 + Av-2) + Pr(N - 1 -> N - l)(N - 1 + Av- 

7TJV-1 

(B.3) 

where A is defined as the expected cumulative weight of the number of descendants 
from a mutant lineage starting with i mutants, conditioned on eventual extinction of 



the lineage (Weissman et al. 2009). This weight represents the mutational opportunity 



for the lineage of primary mutants and is scaled by the population size. The boundary 



1G 



conditions are that D = while does not need to be defined since no transition to 
state N is possible. In the neutral case where r = 1 this reduces to 



A = i(l + D„) + (l - ^^r^) (1 + A) + ^(1 + D 2 ) (B.4) 

A = + A-,) + (i - (« + A) + ^9^> (. + %■) 

Dm . m^V, - 1 + D „_ 2 ) + (l - (N - 1 + Djv_i). 



The solution = iN /2 satisfies system (B.4 1. Because generations are scaled in terms 



of N time steps in the Moran model we can say that on average a single mutant produces 
an effective number of N/2 descendants before going extinct. That is, the sum of the 
time that primary mutant descendants are alive is N/2 generations. 

In the case where r < 1 we can write the system for finite population size but I have 
found no simple way of expressing its solution. In the limit of large population size we 
can make use of the fact that a lineage starting with i mutants must produce i times 
as many descendants as a lineage starting with 1 mutant. I define D* as the number of 
descendants produced scaled to the population size such that D s ; = D*N. Solving for 
D* +1 and taking the limit as ./V — > oo gives 

_ (r + !)£>* -DU-1 
u *+i - ~ ■ 

The additional conditions that = " P '~ 1 = D ) + , 1 implies that 

Thus, for the Moran model, the scaled number of mutants produced by a single mutant 
approaches 1/(1 — r). For finite populations, the value is within 1% of this limit for 
populations larger than 100 when r < 0.75. 

Appendix B.2. Wright- Fisher populations 

For Wright-Fisher populations I calculate the number of descendants when the pri- 
mary mutation is neutral. Consider a population with N haploid genotypes. The first 
step equations are 

N 

D i = ^2Pr(i^j)^-(i + D j ), 
— 

3=0 



Pr(i -> j) = 



N\ ( i V {N-i\ N - j 



j J \N J V N 

N -i 



N 

Because the mutant allele is neutral there is no effect on mean fitness as the number 
of primary mutants changes in the population. This means that, conditioned on non- 
fixation of the primary mutation, the number of descendants left by i primary mutants 
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must just be i times the number left by 1 primary mutant. Thus Di = D\i. Inserting 
this relationship back into the system of equations gives 

which is can be written as the sum of terms involving the first and second (non-central) 
moments of the binomial distribution with parameters N and i/N. Simplification gives 

/ DUN -I) 
Di* = t(l+ V 



N 

Solving for D\ gives D 1 — N. Thus 

Di = Ni. (B.6) 

So for a haploid Wright-Fisher population the mean number of neutral mutants spawned 
by a lineage destined to extinction is equal to the number of haploid genomes present in 
the population. 

For deleterious mutations it is not possible to write the conditional process for finite 
populations because no closed-from solution for the fixation probabilities is available. 
Instead, I calculate the average number of descendants in a large population using the 
Poisson approximation. The first step equations are 

N 

Di = ^2,Pr{i -> + Dj), 

3=0 

e- ir (ir) j 



Pr(i — > j) = 



In very large populations the branching approximation applies and Di = iDi, yielding 

Di = D\i = i + D\ir. 

Solving for D\ shows that 

A - Y~ — ■ (B.7) 
1 — r 
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probability that no successful secondary mutants are spawned from a lineage descending from 
a single primary mutant 
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approximate a derived by Iwasa et al. (2004) 


Tm 


for the Moran model, the probability that a single primary mutant will produce a lineage 
that produces a successful secondary mutant. 


Twf 


for the Wright-Fisher model, the probability that a single primary mutant will produce a lineage 
that produces a successful secondary mutant. 
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Figure B.l: The solution for the probability of tunneling in large Wright-Fisher populations. Panel 
A shows the probability of tunneling as a function of uj and 1 — r (note that th e selection coefficient 
s = 1 — r). When 1 — r ~ the curve is approximately given by equation \21\ . As 1 — r increases 
the probability of tunneling approaches that given by equation | |25| . Panel B shows the decrease of the 
probability of tunneling as 1 — r goes up in blue (ui = 0.00001). The red curve shows the approximate 
value for small y/Zj relative to 1 — r from equation < |25| . For this value of lu, the two curves are within 
5% once 1 - r > 0.01. 
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Figure B.2: The probability of tunneling in multistep pathways under the Moran model can be found 
by recursively applying the formula for the pro babi lity of tunneling. The diagonal line is shown in black 
along with the recursion formula from equation ( |2(i| for three different values of r. In all cases, fi = 10 — 6 . 
The green curve is for r = 1, where each intermediate mutation has the same fitness as the ancestral 
allele. The graph can be cobwebbed as shown by the black arrows. The initial value is the probability 
that the beneficial mutant at the end of the pathway becomes fixed, starting from a single individual. 
Regardless of the initial condition, the probability of tunneling through a long pathway converges on the 
value where the green curve crosses the black line. The blue curve is for slightly deleterious intermediate 
mutants with r = 0.9999995. This still crosses the black line at a positive value. The red curve is for 
a lower value of r = 0.999998. The red curve is always below the black line, so longer pathways have 
essentially no chance of fixing via tunneling. 



22 



true total probability 




minimum maximum 



Number of primary mutants 



Figure B.3: A schematic illustration of the relationship between the probability of tunneling and the 
mean number of primary mutants. The probability that no successful secondary mutants are produced 
is a decreasing function of the primary mutant lineage size with positive second derivative. Assume 
that the primary mutant lineage takes on two possible values with equal probability (the minimum and 
maximum). The red lines map from the primary mutant lineage size and have a mean given by the red 
dashed line. The probability that no successful secondary mutant is produced from events having the 
mean lineage size is shown where the red dashed line meets the vertical axis. The true total probability 
that no successful secondary mutants are produced is found by averaging the probabilities in the two 
different lineage sizes which can be visualized by finding the midpoint between the blue horizontal lines. 
Thus, the true probability that no secondary mutants are produced is higher than that found from the 
mean population size. This means that the probability of tunneling is smaller than would be expected 
based on mean population size alone. 
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Figure B.4: Predicted and observed waiting times for tunneling under the Moran model. The black points 
show the waiting time in units of generations for a successful secondary mutation to arise along with the 
95% confidence intervals. Each parameter set was simulated 10 4 times. The yellow curve represents the 
expected waiting time when the tunneling pathway is ignored. The orange curve represents the Iwasa 
et al. (2004) approximation (their equation 9), while the red curve represents my new approximation. 
The green curve was constructed by numerically solving the system for finite population size. For large 
population size, the two approximations agree. The parameter values were chosen to have extremely 
beneficial secondary mutants so that the tunneling effect is large, fi\ = 10~ 5 , ^2 = 10 — 2 , r = 1, and 
a = 100. The simulations were stopped when the advantageous secondary mutation became fixed or 
there were more than 200 secondary mutants (the probability of loss after that point is approximately 
10 — 400 ). Similar results are obtained for r < 1 (not shown). 
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Figure B.5: Predicted and observed waiting times for tunneling under the Wright-Fisher model. The 
black points show the waiting time for a successful secondary mutation to arise along with the 95% 
confidence intervals. Because the Wright-Fisher simulations take much more time than the Moran 
model simulations, each parameter set was simulated only 1000 times, except for TV = 10 7 which was 
only simulated 100 times. The yellow curve represents the expected waiting time when the tunneling 
pathway is ignored. The orange curve represents the formula used by Lynch and Abegg based on the 
Iwasa et al. (2004) approximation, while the red curve is based on my Wright-Fisher approximation. 
The green curve was constructed by numerically solving the system for finite population size. Because 
this requires solving a non-sparse matrix equation that is the size of the population, I was only able to 
numerically solve the finite population size model for population size less than 10 4 . Panel B shows the 
graph for these population sizes in more detail. For small population sizes, the simulated data match the 
numerical prediction quite well. For population sizes larger than about 10 4 the observed values agree 
with my Wright-Fisher approximation. Parameter values are fi\ = 10 — 8 , fi2 = 10 -5 , r = 1, a = 1.01. 
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